1/1 


MICROCOPY  resolution  test  chart 

_  frATIQftAl  EkJRLAO  Of  STANDARDS  J9%3  A 


Special  Report  87-8 

June  1987 


US  Army  Corps 
of  Engineers 


OTIC  FILE  CO'Pl 


Cold  Regions  Research  & 
Engineering  Laboratory 


Tracking  two-dimensional  freezing  front 
movement  using  the  complex  variable 
boundary  element  method 


Ted  Hromadka 


DTIC 


ELECTE 

k  AUG  1  4  1987 


T 


Approved  for  public  release;  distribution  Is  unlimited. 


87  a  1 


A  sr  o 


L*  ’if. 


it,  |1 


|l. 


at 


'•'■  'V 


CURITY  CLASSIFICATION  OF  THIS  PAGE 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No  0704-0188 
Exp  Oare  lun  )0.  1986 


la.  REPORT  SECURITY  CLASSIFICATION 

Unclassified 


lb  RESTRICTIVE  MARKINGS 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


2b.  DECLASSIFICATION  /  DOWNGRADING  SCHEDULE 


3.  DISTRIBUTION  /  AVAILA8ILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlmited. 


4.  PERFORMING  ORGANIZATION  REPORT  NUMSER(S) 


S  MONITORING  ORGANIZATION  REPORT  NUM8ER(S) 

Special  Report  87-8 


6a  NAME  OF  PERFORMING  ORGANIZATION  1 6b.  OFFICE  SYM80L  I  7a  NAME  OF  MONITORING  ORGANIZATION 


Williamson  &  Schmid 


(If  applicable) 


6c.  ADDRESS  (City.  State,  and  ZIP  Code) 

Irvine,  California  92714 


U.S.  Army  Cold  Regions  Research 
and  Engineering  Laboratory 


7b.  ADDRESS  (City.  State,  and  ZIP  Code) 


Hanover,  New  Hampshire  03755-1290 


8a.  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 


8b  OFFICE  SYMBOL 
(If  applicable) 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


P.O.  84-M-1691 


8c.  ADDRESS  (City.  State,  and  ZIP  Code) 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

TASK 

NO 

NO 

WORK  UNIT 
ACCESSION  NO 


CW1S  31711 


11.  TITLE  (Include  Security  Classification) 

Tracking  Two-dimensional  Freezing  Front  Movement  Using  the  Complex  Variable 
Boundary  Element  Method _ 

12.  PERSONAL  AUTHOR(S) 

Ted  Hrodmadka _ _ _  _ _ _ _ 

13a.  TYPE  OF  REPORT  3b  TIME  COVERED  |1 4  DATE  OF  REPORT  (Year.  Month,  Day)  |l 5  PAGE  COUNT 

from _ to _  June  1987  64 


16.  SUPPLEMENTARY  NOTATION 


Ue  C  V'Bft'V 


17 

COSATI  COOES  | 

FIELD 

GROUP 

SUB-GROUP 

18.  SUBJECT  TERMS  ( Continue  on  reverse  if  necessary  and  identify  by  block  number) 

Computer  models  Soils 

Frozen  soils  \  Soil  water  ^ _ _ 

Heat  flux  \  Water  .y'f'kt  tVftF P  i ) 


mating  the  location  of  the  freezing  front  in  soil- water  phase  change  problems.  ^Fbis^computer  program<^_ 
®VBFR4>  is  based  on  the  following  major  assumptions:  1)  the  problem  is  twopdimensional;  2)  the  entire  soil 
system  is  homogeneous  and  isotropic;  3)  the  problem  thermal  boundary  conditions  are  constant  values  of 
temperature  (or  stream  function);  4)  soil/Water  flow  effects  are  neglected  (the  problem  is  strictly  geother¬ 
mal);  5)  all  heat  flow  from  the  freezing  front  is  within  the  control  volume,  there  is  no  heat  flux  associated 
with  the  freezing  front  from  exterior  of  the  control  volume;  and  6)  the  freezing  front  movement  is  slow 
enough  that  heat  flux  along  the  moving  boundary  can  be  determined  by  assuming  steady  state  heat  flow 
conditions  for  small  durations  of  time  (i.e.,  timesteps).  The  CVBEM  is  used  to  model  the  thermal  regime 
of  the  soil  system.  The  theory  a»  i  development  of  the  CVBEM  are  given  in  CRREL  Internal  Xeport  969,  — ■ 
"Complex  Variable  Boundary  Elements  in  Engineering,"  by  Hromadka.  Because  the  numerical  technique  is 
a  boundary  integral  approach,  the  control  volume  thermal  regime  is  modeled  with  respect  to  the  boundary 
values,  and,  therefore,  the  CVBFR1  data  entry  requirements  are  significantly  less  than  those  usually  re- 
quired  of  domain  methods.such  as  finite-differences  orfinlte-elements>  Soil-water  phase  change  along  the 

20  DISTRIBUTION  /  AVAILA8ILITY  OF  ABSTRACT  21  ABSTRACT  SECURITY  CLASSIFICATION 

EDGnclassified/unlimited  □  same  as  rpt  □  one  users  Unclassified  _ _ 


20  DISTRIBUTION /AVAILABILITY  OF  ABSTRACT  21  ABSTRA 

EDOnclassified/unlimited  □  same  as  rpt  □  one  users  Unclai 

22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL  22b  TElEPH 

__££3-6 

DD  FORM  1473,  84  MAR  83  APR  edition  may  be  used  until  exhausted 

All  other  editions  are  obsolete 


22b  TELEPHONE  (Include  Area  Code)  I  22c  OFFICE  SYMBOL 


603-646-4335 


CECRL-EE 


SECURITY  CLASSIFICATION  OF  THIS  PAGE 

Unclassified 


.V 


■V/.V  ■ 


rv. 


* 

# 


Unclassified 


19.  Abstract  (cont'd) 

X  K 

-^freezing  front  is  modeled  as  a  simple  balance  between  computed  heat  flux  and  Wte  evolution  of 
soil-water  volumetric  latent  heat  of  fusion.  To  model‘£hff,displacement  of  the  freezing  front, 
program  CVBFR1  provides  two  options:  4f*aisplace  the  freezing  front  coordinates  with  respect 
to  changes  in  the  y-coordinate  only,^2^dispiace-the-freezing  front  coord inate^nvith  respect  to 
a  vector  normal  to  the  freezing  front  boundary. 


Unclassified 


PREFACE 

This  report  was  prepared  by  Ted  Hromadka,  Director  of  Water  Re¬ 
sources,  Williamson  and  Schmid.  The  work  was  performed  for  CRREL 
under  Contract  84-M-1691  and  was  funded  by  the  Directorate  of  Civil 
Works,  Office  of  the  Chief  of  Engineers,  under  Civil  Works  Order  No. 
CWIS  31711,  Time  Rate  and  Magnitude  of  Degradation  of  Permafrost. 

Critical  reviews  of  the  report  were  furnished  by  Dr.  Richard 
Berg  and  Francis  Sayles,  the  project  monitor. 

The  contents  of  this  report  are  not  to  be  used  for  advertising  or 
promotional  purposes.  Citation  of  brand  names  does  not  constitute 
an  official  endorsement  or  approval  of  the  use  of  such  commercial 
products. 


Accesston  ?or 

HTIS  GV.fcl 
DTIC  T AS 
Unannounced  □ 

Justification _ 


By - 

Distribution/ 

Availability  Codes 
[Avail  and/or 
j  Special 


iDist 


TABLE  OF  CONTENTS 


INTRODUCTION 

Overview 

Objectives  of  Report 
Report  Organization 
Report  Preparation 

MODELING  APPROACH 

Introduction 

Heat  Flow  Model 

Phase  Change  Model 

Program  CVBFR1  Characteristics 

PROGRAM  CVBFR1 
Introduction 
Problem  Set-Up 
Input  Data 
Application 

APPENDIX 

COMPLEX  VARIABLE  BOUNDARY  ELEMENT  METHOD 

THE  APPROXIMATIVE  BOUNDARY 
FOR  CVBEM  ERROR  ANALYSIS 

PROGRAM  CVBFR1  FORTRAN  LISTING 


LIST  OF  FIGURES 


Typical  Two- Phase  Problem  Definition 
Typical  Roadway  Embankment  Problem 
Nodal  Point  Placement  and  Boundary  Conditions  for 
Fig.  2.2  Problem 

Normal  Vector  Coordinate  Displacement  Model 
Program  CVBFR1  Boundary  Condition  Characteristics 
CVBFR1  Modeling  Procedure 

Example  Problem  Roadway  Embankment  Discretized  into 
Finite  Elements 

Nodal  Point  Numbering  for  4  CVBEM  Nodal  Densities 
Comparison  of  CVBEM  Model  Results  in  Predicting  Freezing 
Front  Location 

Initial  Conditions  and  Cross-Section  Locations 
Comparison  of  CVBEM  and  NDI  Modeling  Results 

Simply  Connected  Domain  fi  with  Simple  Closed 
Contour  Boundary  T 

r  Discretized  into  m  Boundary  Elements 

(k  +  1)  -  Node  Boundary  Element  Ti  Nodal  Definitions 

Branch  Cut  of  LN(Z-  $")  Function 

The  Analytic  Continuation  of  u5  ( Z  )  to  the 
Exterior  of  QuT 

Application  Problem  Geometries  and  Solutions  for 
Tern  perature,  ( x ,  y  ) 

Approximative  Boundaries  for  Three  Nodal 
Point  Distributions 

Approximative  Boundaries  for  the  Five  Nodal 
Point  Distributions 


.  /.  A  V. 


1.  INTRODUCTION 


1.0  OVERVIEW 


The  Complex  Variable  Boundary  Element  Method  or  CVBEM  is  used  to  develop 
a  computer  model  for  estimating  the  location  of  the  freezing  front  in  soil- 
water  phase  change  problems.  This  computer  program,  CVBFR1 ,  is  based  on 
the  following  major  assumptions: 

(i)  the  problem  is  two-dimensional, 

(ii)  the  entire  soil  system  is  homogeneous  and  isotropic, 

(iii)  the  problem  thermal  boundary  conditions  are  constant  values 
of  temperature  (or  stream  function), 

( i v )  soil-water  flow  effects  are  neglected  (the  problem  is  strictly 
geothermal ) , 

(v)  all  heat  flow  from  the  freezing  front  is  within  the  control 
volume;  there  is  no  heat  flux  associated  with  the  freezing  front 
from  exterior  of  the  control  volume, 

(vi)  the  freezing  front  movement  is  sufficiently  slow  such  that 
heat  flux  along  the  moving  boundary  can  be  determined  by 
assuming  steady  state  heat  flow  conditions  for  small  durations 
of  time  (i.e.,  timesteps). 

The  CVBEM  is  used  to  model  the  thermal  regime  of  the  soil  system. 

The  theory  and  development  of  the  CVBEM  is  given  in  Hromadka  (1984,  1987). 
Because  the  numerical  technique  is  a  boundary  integral  approach,  the 
control  volume  thermal  regime  is  modeled  with  respect  to  the  boundary 
values  and,  therefore,  the  CVBFR1  data  entry  requirements  are  signifi¬ 
cantly  less  than  that  usually  required  of  domain  methods  such  as  finite- 
differences  or  finite-elements. 
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Soil-water  phase  change  along  the  freezing  front  is  modeled  as  a 
simple  balance  between  computed  heat  flux  and  the  evolution  of  soil- 
water  volumetric  latent  heat  of  fusion.  To  model  the  displacement  of 
the  freezing  front,  program  CVBFR1  pro’ ides  two  options: 

(i)  displace  the  freezing  front  coordinates  with  respect  to 
changes  in  the  y-coordinate  only, 

(ii)  displace  the  freezing  front  coordinates  with  respect  to 
a  vector  normal  to  the  freezing  front  boundary. 

1.1  OBJECTIVES  OF  REPORT 

The  objectives  of  this  report  are  threefold: 

(1)  Provide  background  information  regarding  the  CVBEM  and  the 
soil -water  phase  change  model  used  in  program  CVBFR1. 

(2)  Provide  documentation  for  the  data  entry  sequence  assoc¬ 
iated  with  program  CVBFR1 . 

(3)  Because  the  CVBEM  results  in  a  small  FORTRAN  computer 
programming  effort,  provide  the  CVBFR1  computer  code  as 
an  appendix  to  this  report. 

1.2  REPORT  ORGANIZATION 

This  report  is  organized  into  4  chapters  and  3  appendices  as  follows: 


SECTION 
Chapter  1 
Chapter  2 


Chapter  3 
Appendix  A 
Appendix  B 


Appendix  C 


Introduction 

Modeling  approach.  Presents  heat  flow  model 
(CVBEM)  and  phase  change  approximation. 

Data  input  requirements  for  program  CVBFRl . 

Background  development  of  the  CVBEM. 

Background  development  of  the  approximative 
boundary  technique  to  evaluate  CVBEM  approxi¬ 
mation  error. 

Program  CVBFRl  source  code. 
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2.  MODELING  APPROACH 


2.0  INTRODUCTION 

The  use  of  the  Complex  Variable  Boundary  Element  Method  to  model 
soil-water  phase  change  effects  is  a  new  numerical  approach  to  this  class 
of  problems.  In  previous  work,  Hromadka  and  Guymon  (1982)  applied  the 
complex  variable  boundary  element  method  (CVBEM)  to  the  problem  of  predict¬ 
ing  freezing  fronts  in  two-dimensional  soil  systems.  Hromadka  et  al . 

(1983)  subsequently  compare  the  CVBEM  solution  to  a  domain  solution  method 
and  prototype  data  for  the  Deadhorse  Airport  runway  at  Prudhoe  Bay,  Alaska. 

In  another  work,  the  model  is  further  extended  to  include  an  approximation 
of  soil-water  flow  (Hromadka  and  Guymon,  1984a).  In  contrast  to  the  CVBEM 
approach,  an  example  in  the  use  of  real  variable  boundary  element  methods 
(Brebbia,  1978)  in  the  approximation  of  such  moving  boundary  phase  change 
problems  and  a  review  of  the  pertinent  literature  is  given  in  O'Niell  (1983). 

Hromadka  and  Guymon  (1984b)  develop  a  relative  error  estimation  scheme 
which  exactly  evaluates  the  relative  error  distribution  on  the  problem 
boundary  that  results  from  the  CVBEM  approximator  matching  the  known  boun¬ 
dary  conditions.  This  relative  error  determination  is  used  to  add  or  delete 
boundary  nodes  to  improve  accuracy.  Thus,  the  CVBEM  permits  a  direct  and 
immediate  determination  of  the  approximation  error  involved  in  solution  of 
an  assumed  Laplacian  system.  The  modeling  accuracy  is  evaluated  by  the 
model-user  in  the  determination  of  an  approximative  boundary  upon  which 
the  CVBEM  provides  an  exact  solution.  Although  inhomogeneity  (and  aniso¬ 
tropy)  can  be  included  in  the  CVBEM  model,  the  resulting  ful ly-populated 
matrix  system  quickly  becomes  large.  Therefore  in  this  work,  the  domain 


is  assumed  homogeneous  and  isotropic  except  for  differences  in  frozen  and 
thawed  conduction  parameters  for  freezing  and  thawing  problems,  respectively. 

A  major  benefit  in  the  use  of  the  CVBEM  over  other  numerical  methods 
(including  real  variable  boundary  element  methods  and  domain  methods  such 
as  finite-differences  and  finite-elements)  is  the  accurate  and  easy-to-use 
"approximative  boundary"  error  evaluation  technique.  Other  numerical  methods 
can  be  evaluated  for  modeling  error  (where  exact  mathematical  solutions 
do  not  exist)  by  increasing  nodal  point  densities  and  comparing  the  result¬ 
ing  changes  in  predicted  nodal  values  of  the  governing  equation's  state 
variable.  In  contrast,  the  CVBEM  approximative  boundary  error  evaluation 
technique  is  simply  the  process  of  locating  the  (x,y)  points  where  the 
CVBEM  approximation  function  meets  the  specified  boundary  condition  values 
(the  approximative  boundary),  and  comparing  the  resulting  plot  to  the  true 
problem  boundary. 

A  major  benefit  for  using  the  CV8EM  error  evaluation  technique  is 
that  highly  accurate  solutions  for  two-dimensional  potential  problems  can 
be  obtained.  Often,  the  CVBEM  approximation  analysis  is  terminated  when 
the  approximative  boundary  differs  from  the  true  problem  boundary  to 
within  the  construction  tolerance  of  the  project,  resulting  in  an  exact 
CVBEM  model  of  a  probable  constructed  version  of  the  engineered  plan  draw¬ 
ings.  Consequently  the  CVBEM  approach  can  be  used  directly  in  engineering 
applications,  or  used  to  provide  a  wide  range  of  highly  accurate  approxi¬ 
mations  for  two-dimensional  phase  change  problems  (where  the  freezing 
front  movement  is  slow;  see  section  2.2)  for  checking  modeling  results 
produced  by  other  numerical  methods. 


2.1  HEAT  FLOW  MODEL 


For  a  wide  range  of  soil  freezing  (or  thawing)  problems,  the  freezing 
front  movement  is  sufficiently  slow  such  that  the  governing  heat  flow  equation 
can  be  modeled  using  a  timestepped  steady  state  heat  flow  approximation. 

That  is  for  small  durations  of  time,  the  heat  flux  along  the  freezing  front 
can  be  computed  assuming  the  temperature  distribution  within  the  frozen 
(or  thawed)  regions  are  potential  functions  (i.e.,  the  Laplace  equation 
applies).  Figure  2.1  illustrates  a  typical  two-phase  problem  definition 
where  the  heat  flow  model  solves  for  heat  flux  along  the  freezing  front  by 
solving  the  Laplace  equation  (by  use  of  potential  functions)  in  both  the 
frozen  and  thawed  regions. 

To  develop  mathematical  models  of  the  Laplace  equation  in  each  region, 
a  CVBEM  approximator  is  generated  which  matches  specified  boundary  conditions 
of  either  temperature  or  flux  at  nodal  point  locations  on  the  problem 
boundary  and  freezing  front.  The  CVBEM  approximator  exactly  satisfies  the 
Laplace  equation;  consequently  there  is  no  modeling  error  in  solving  the 
governing  Laplace  equation  (heat  flow  model),  there  is  only  error  in  matching 
the  boundary  conditions  continuously.  Figure  2.2  shows  an  example  roadway 
problem  where  the  freezing  front  is  initially  located  some  known  distance 
below  the  surface.  Boundary  conditions  for  the  example  problem  and  a  nodal 
point  placement  scheme  are  shown  in  Fig.  2.3. 

The  heat  flow  model  in  CVBFR1  develops  a  CVBEM  potential  function  which 
satisfies  the  Laplace  equation  within  the  boundary  of  Fig.  2.3.  Appendix  A 
provides  a  brief  review  of  the  CVBEM  numerical  approach,  and  Appendix  B 
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FIG.  2.1  TYPICAL  TWO-PHASE  PROBLEM  DEFINITION 
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FIG.  2 .2  TYPICAL  ROADWAY  EMBANKMENT  PROBLEM 
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provides  a  review  of  the  approximative  boundary  error  evaluation  technique 
used  to  develop  more  accurate  CVBEM  approximations.  The  usual  modeling 
procedure  is  to  use  the  approximative  boundary  technique  to  analyze  the 
initial  condition  CVBEM  model.  After  the  analyst  is  satisfied  with  the 
CVBEM  approximator  and  its  associated  level  of  accuracy  then  the  CVBFR1 
program  is  executed  to  model  the  freezing  front  evolution. 

2.2  PHASE  CHANGE  MODEL 

For  each  timestep,  a  CVBEM  approximator  is  generated  by  program  CVBFR1 
based  on  the  problem  geometry  and  boundary  conditions.  Heat  flux  is  computed 
along  the  freezing  front  using  the  CVBEM  approximation  stream  function 
values.  The  heat  flux  estimates  are  assumed  to  directly  equate  to  the 
rate  of  freezing  (or  thawing)  of  a  volume  of  soil  at  the  freezing  front. 
Consequently,  a  freezing  process  for  the  example  of  Fig.  2.3  results  in  a 
downward  migration  of  the  freezing  front  such  that  the  product  of  the  time- 
step  and  heat  flux  equals  the  latent  heat  evolved  by  the  change  in  freezing 
front  coordinates. 

Two  freezing  front  displacement  models  are  available  in  program  CVBFR1: 

(1)  All  displacement  occurs  in  the  vertical  direction.  This 
simplified  model  is  generally  appropriate  for  many  road¬ 
way  problems. 

(2)  Displacement  computed  based  on  an  outward  normal  vector. 

This  model  is  the  most  accurate,  but  requires  more  compu¬ 
tational  effort  than  the  vertical  displacement  model. 

Figure  2.4  shows  the  nodal  point  displacement  in  a 
direction  which  balances  the  angles  to  go  between  the 
normal  vector  and  boundary  elements. 
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FIG.  2.3  NODAL  POINT  PLACEMENT  AND  BOUNDARY  CONDITIONS 
FOR  FIG.  2.2  PROBLEM 
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FIG.  2.4  NORMAL  VECTOR  COORDINATE  DISPLACEMENT  MODEL 
(not*  balanced  angles  for  oach  normal  vector) 


2.3  PROGRAM  CVBFR1  CHARACTERISTICS 


Class  of  Problems  Modeled 

Program  CVBFR1  may  be  used  to  model  soil-water  freezing  (or  thawing)  in 
two-dimensional,  homogeneous,  isotropic  domains.  As  illustrated  by  the 
example  problem  of  Figs.  2.2  and  2.3,  only  one  region  is  modeled  ( i . e . , 
either  entirely  frozen  or  entirely  thawed)  and  the  freezing  front  forms  part 
of  the  control  volume's  boundary.  For  example,  program  CVBFR1  may  be  used  to 
study  the  freezing  front  advancement  into  a  soil  system  where  the  soil  sys¬ 
tem  is  initially  close  to  the  freezing  point  depression  temperature,  and 
negligible  heat  flow  to  the  freezing  front  is  contributed  from  the  under¬ 
lying  soil  system.  A  schematic  of  the  problem  domain  and  boundary  conditions 
used  in  CVBFR1  are  illustrated  in  Fig.  2.5.  Another  characteristic  of 
CVBFR1  is  that  the  boundary  conditions  of  the  problem  are  held  constant 
for  the  entire  simulation.  Additionally,  the  initial  conditions  of  the  prob¬ 
lem  are  assumed  to  be  near  steady  state  with  the  freezing  front  specified 
some  distance  below  the  top  of  the  control  volume  boundary  (control 
surface) . 

The  CVBFR1  Modeling  Procedure 

The  modeling  procedure  used  in  the  CVBFR1  program  is  shown  schematically 
in  Fig.  2.6  for  the  case  of  a  soil  freezing  problem.  It  is  assumed  in 
Fig.  2.6  that  the  analyst  has  developed  a  good  CVBEM  approximator  for  the 
initial  conditions  of  the  problem  by  using  the  approximative  boundary 
technique  (Appendix  B)  to  locate  nodal  points  on  the  problem  boundary. 
Typically,  the  most  difficult  modeling  problem  occurs  when  the  freezing  front 
is  closest  to  the  top  of  the  problem  boundary  such  as  shown  in  Fig.  2.3. 
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DEVELOP  A  CVBEM  APPROXIMATOR  BASED 
ON  BOUNDARY  COORDINATES  AND  BOUN¬ 
DARY  CONDITIONS 


T=  -  IO°C 


CALCULATE  HEAT  FLUX  VALUES  ALONG 
THE  FREEZING  FRONT 


DISPLACE  NODAL  COORDINATES  ALONG 
FREEZING  FRONT  BASED  ON  HEAT  EVOLVED, 
AND  VOLUMETRIC  LATENT  HEAT  OF  FUSION 
FOR  SOIL-WATER  MIXTURE 


FIG.  2.6  CVBFR1  MODELING  PROCEDURE 


12 


Consequently,  the  CVBEM  nodal  placement  should  be  concluded  based  on  the 
problem's  smallest  anticipated  distance  to  the  freezing  front.  For  example, 
the  roadway  problem  shown  in  Fig.  2.3  spans  a  width  of  50  m;  the  corre¬ 
sponding  distance  to  the  freezing  front  for  initial  conditions  (freezing 
problem)  is  assumed  to  be  0.25m. 


3.  PROGRAM  CVBFR1 


3.0  INTRODUCTION 

CVBFR1  is  a  CVBEM  program  with  the  capability  of  estimating  the  moving 
position  of  a  slow-moving  freezing  front  in  soils.  The  CVBFR1  program  uses 
either  subroutine  FRT1  or  FRT2  to  estimate  the  displacement  of  the  freezing 
front  where  subroutine  FRT1  is  based  upon  a  vertical  shifting  and  FRT2  uses 
the  outer  normal  direction  to  calculate  the  change  in  nodal  point  coordinates. 

3.1  PROBLEM  SET-UP 

The  problem  domain  is  assumed  to  be  a  homogeneous  isotropic  soil  mix¬ 
ture  enclosed  by  the  problem  boundary.  Nodal  points  are  located  on  the  prob¬ 
lem  boundary  and  are  numbered  in  sequence  in  a  counterclockwise  direction 
from  1  to  NNOD. 

Nodal  points  are  generally  placed  closer  together  near  angle  points  of 
the  problem  boundary,  or  where  boundary  condition  values  (or  types  of 
boundary  conditions)  change.  This  increase  in  nodal  density  reduces  the 
error  in  integrating  a  trial  function  (straight  line  interpolation  functions 
are  used  in  CVBFR1)  which  becomes  inaccurate  near  singularities  of  the  po¬ 
tential  function,  temperature. 

The  product  of  the  latent  heat  of  fusion  for  soil -water  and  the  uniform 
soil  porosity  value  is  used  as  the  volumetric  latent  heat  of  fusion  for  the 
soil-water  (or  soil-ice)  mixture.  The  thermal  conductivity  value  is  used 
to  estimate  the  normal  heat  flux  values  along  the  freezing  front. 


3.2  INPUT  DATA 


Input  data  for  program  CVBFR1  is  as  follows: 

VARIABLE  DATA  FILE  LINE 

KODE  Line  1 

NNOD,  NFRS,  NFRE  Line  2 

COND,  XLAT,  POR  Line  3 

DELT,  SIMUL,  OUT,  ID  Line  4 

X(I),  Y( I ) ,  KTYPE(I),  VALUE(I);  1=1  to  NNOD  Line  5 


X(NNOD) ,  Y(NNOD) ,  KTYPE(NNOD) ,  VALUE (NNOD);  Line  NNOD  +  4 
(END  OF  FILE) 
where: 


VARIABLE 

KODE 


NNOD 

NFRS 

NFRE 

COND 

XLAT 

POR 

DELT 

SIMUL 


1,  For  vertical  displacement  of  freezing 
front  coordinates 

2,  Use  outward  normal  vector  to  estimate 
nodal  point  displacements 

Total  number  of  nodes  on  boundary 

First  node  number  of  the  freezing 
front  contour 

Last  node  number  of  the  freezing 
front  contour 

Thermal  conductivity  of  a  homogeneous 
isotropic  soil  mixture 

Latent  heat  of  fusion  for  soil -water 

Porosity  of  soil 

Increment  for  time  advancement  model 
Total  simulation  time 


OUT 


Output  period 


ID  =  0,  Detailed  output  (see  Example  1) 

1,  Summary  output  (see  Example) 

X(I),  Y(I)  =  (x,y)  coordinates  of  node  I  in  first 

quadrant 

KTYPE(I)  =  1,  Prescribed  temperature  value 

2,  Prescribed  stream  function  value 

3,  Prescribed  flux  value 

VALUE(I)  =  Prescribed  value  according  to 

KTYPE(I).  For  efflux,  VALUE(I)  = 

efflux/conductivity 

Note:  The  units  of  XLAT,  COND,  DELT,  SIMUL,  and  OUT  should  be 
consistent. 


3.3  APPLICATION 


Example  1:  Computing  the  Freezing  Front  Location  in  a  Roadway  Embankment 
A  roadway  embankment  (Fig.  3.1)  problem  is  used  to  illustrate  the 
application  of  program  CVBFR1.  The  input  data  and  program  output  (in 
English  units)  for  the  example  problem  is  provided  in  the  following: 
(note  that  the  first  line  is  a  "1"  or  "2"  for  using  subroutines  FRT1  and 
FRT2,  respectively): 


Fig.  3.1b  Nodal  Point  Numbering  fo 


«•>  «I<|A 


,*  i.ii  M  M  lJ 


PROGRAM  CVBFR 1  Data  Input 
(Example  Problem) 


KODE  =  1  or  2  (see  text) 
JO  1  15 
10  30  .4 
6  120  120  0 
0  10  2  o 
25  LO  1  O 
4  9  10  1  0 
49.9  io  l  0 

50  10  1  O 

50.1  9.95  1  0 

51  9.5  1  0 
*0  5  l  0 

59  .510 
.*>9.9  .05  1  0 
70  0  1  0 

70.1  0  I  0 
VI  0  1  0 
95  0  1  0 
120  0  1  0 
120  1  1  -10 
95  1  1  -10 
/l  1  1  -10 

70.1  1  1  -10 
70  1.  1  -10 
<*>9.9  1.05  1  -10 
.*  >9  1.5  l  -10 

60  6  1  -10 

51  10.5  1  -10 

50.1  10.95  1  -10 
50  l l  1  -10 
49.9  11  1  -10 

49  l l  1  -10 


25  11  i  -10 
0  11  1  -IO 
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The  computer  modeling  results  using  FRT1  (KODE  =  1)  are  as  follows: 


TIME  INCREMENT  =  6.0000 

TOTAL  SIMULATION  TIME  =  120.0000 
CONDUCTIVITY  =  10.0000 

LATENT  HEAT  =  SO . 0000 
POROSITY  =  0.4000 


NODE 

X(I1 

Y  i  I ) 

KT  YPE (  I  > 

VALUE 

ANGLE ( I ) 

NO . 

1  =S V  • 2  =  SF 

3-EFFLUX 

1 

0. 00000 

1 o , ooooo 

2 

0 . ooooo 

90.00 

25 . Ooooo 

lo . ooooo 

1 

0 , ooooo 

180.00 

3 

49 . OoOOO 

lo. ooooo 

1 

0 . ooooo 

180.00 

4 

49 . voooo 

10 . ooooo 

1 

0 , ooooo 

180.00 

5 

50 . (>0000 

l 0 . OOoOo 

1 

0 . ooooo 

206.57 

6 

50. 10000 

9.95000 

1 

0 . 0O0OO 

180 . 00 

7 

51 .00000 

9 . 50000 

1 

O . OOOOO 

180 .  (id 

8 

60.00000 

5  •  OOoOO 

1 

0 . ooooo 

180.00 

9 

69 . OOOOO 

0 . 50000 

1 

0 . ooooo 

180.00 

10 

69 . VOOOO 

0 , 05000 

1 

()  .  ooooo 

130.00 

1  l 

70 , OOOOO 

0 , ooooo 

1 

•  > . ooooo 

153.43 

12 

70  ,  J.OOOO 

0  .  ooooo 

1 

0 . ooooo 

180.00 

13 

71 . ooooo 

0 , ooooo 

1 

0  ,  (>0000 

180.00 

14 

95 . OOOOO 

0 . ooooo 

1 

0 . ooooo 

180.00 

15 

120. OOOOO 

0 . ooooo 

1 

0 . OO000 

90 . 00 

16 

120.00000 

1 . ooooo 

1 

-10.00000 

90 . 00 

17 

95 . OOOOO 

1 . 00OOO 

1 

-10.00000 

1  SO .  00 

18 

7 l . OoOOO 

1 , ooooo 

1 

-  10 . Ooooo 

1 9<> .  00 

19 

70 , lOOOO 

1 . ooooo 

1 

- 10 • ooooo 

180.00 

20 

70 . OOOOO 

l - ooooo 

1 

- io . ooooo 

206.57 

21 

69 . 90000 

1  -  >>5000 

1 

- 10 , oOOOO 

180.00 

69 . OOOOO 

1 .50000 

1 

-10. ooooo 

ISO .  00 

23 

6<> .  OOOOO 

6.00000 

1 

-In. OOOOo 

1 80 . 00 

24 

51.00000 

10. 50000 

1 

-10. ooooo 

180.00 

25 

50.10000 

10.95000 

1 

-10 . 0000O 

180.00 

26 

50.00000 

1 1  ,  ooooo 

1 

-10. ooooo 

153.43 

27 

49.90000 

11. ooooo 

1 

-10. ooooo 

180.00 

28 

49 . OOOOO 

1 1 . ooooo 

1 

-10. ooooo 

180 . 00 

29 

25 . ooooo 

11 . ooooo 

1 

-10 . oOOOO 

180.00 

30 

0 , OOOOO 

11 . ooooo 

1 

- lo . ooooo 

90 . 00 

Cauchy  Program  Results 


TIME  = 

120.0000 

NODE 

STATE 

STREAM 

NUMBER 

VARIABLE 

function 

1 

-0.0645 

0.0000 

2 

o , oooo 

186.8884 

3 

0.0000 

366-0428 

4 

0 . OOOO 

373.0393 

5 

0 . oooo 

373.7330 

6 

o . oooo 

374.6362 

7 

<>.0000 

382.6803 

8 

0.0000 

462.581 1 

9 

o . oooo 

542.9409 

10 

0 . oooo 

550.4897 

11 

0  .  oooo 

551.2713 

12 

o.oooo 

552.0007 

13 

0 , oooo 

558.4281 

14 

0.0000 

739,2386 

15 

o . oooo 

927. 1421 

16 

-JLO.O000 

927.1481 

17 

-lo.oooo 

739.2394 

18 

-10.0000 

558,3662 

19 

-10. oooo 

550.7808 

20 

-10,0000 

549.4852 

21 

-10.0000 

548,0763 

22 

-10.0000 

539.4200 

23 

-10 . oooo 

459-0154 

24 

-10 , oooo 

379.2508 

25 

-10.0000 

372,4201 

26 

-10.0000 

372.0577 

27 

-lo.oooo 

371.7505 

28 

-10.0000 

365.9584 

29 

— 1 O . oooo 

186.8909 

30 

-  to . oooo 

(> .  0024 

CVBEM  Approximation  Function  Nodal  Values: 


NODE 

NUMBER 

1 

2 

3 

A 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 


state 

variable 

-0.0582 
-0.0123 
-0.0574 
-0.1052 
0.0292 
-O. 1064 
-0.0448 
-0.0065 
0.0311 
0.0745 
-o -  0335 
0.0556 
o .  ('>402 
-0.0026 
o.ooi 3 
-  10,0028 
-1 0.021 6 
-9.9618 
-9.9189 
-9 . 8962 
-9.9198 
-9. 9632 
-9.9878 
- i 0  *  0409 
-  1 0.0944 

-10.0653 

-i(l,  0806 
-10. 05 1 0 
-9,9580 
-9.9981 


STREAM 
FUNCTION 
-0.0034 
1 86.8887 
366. 1034 
373. 11 09 
373.7366 
374.5667 
382  <  6476 
462.5815 
542.9185 
550.4381 
551 .2729 
552.0399 
558.4709 
739.2393 
927. 1403 
927.1437 
739.2404 
558.3270 
550.7245 
549.4889 
548. 1329 
539.4475 
459.0158 
379.2812 
372.4865 
372.0609 
371.6945 
365.9055 
186.8916 
-0.0011 
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Nodal  Point  Relative  Error  Values: 

-0.0063  0,0034 

*''•<'>123  -0.0004 

<>.0574  -0.O6O5 

0.1052  -0.0717 

-0.0292  -0.0036 

0.1064  0.0695 

0.0448  0.0327 

0.0065  -0.0004 

-0.0311  0.0225 

-0.0745  0,0516 

0,0335  -0.0016 

-0.0556  -0,0392 

-0.0402  -0.0428 

0.0026  -0.0007 

-0.0013  0.00I8 

0.0028  0.0045 

0.0216  -0.0010 

-0.0382  0.0392 

-0.0811  0,0563 

“0-1038  -0.0037 

-0.0802  -0.0566 

-0.0368  -0.0275 

-o .0122  -0,0005 

0.0409  -0.0304 

0-0944  -0.0664 

0.0653  -0,0031 

0.0806  0.0560 

0 . 0510  0.0529 

“0-0420  -0.0007 

-0.0019  0,0035 

New  Coordinates  of  the  Freezing  Front 


Node 

X-Coord. 

Y-Coord. 

1 

0 , 0000 

9,6542 

6 

24.9999 

9.6543 

3 

48.9923 

9.6363 

4 

49.8163 

9.5782 

5 

49,8736 

9.5383 

6 

49-9673 

9.5380 

7 

50.8326 

9, 1534 

8 

59.832? 

4.6640 

9 

68.8452 

0.1761 

IO 

69.8043 

-O, 2259 

11 

69,9464 

-0.2389 

12 

70,0668 

-0.2815 

13 

70.9948 

-0.3342 

14 

95.0000 

-0,3469 

15 

120  <  OOOO 

-0.3468 

The  output  (summary)  data  using  FRT2  (KODE  =  2)  consists  of: 


TIME  INCREMENT  =  6.0000 

TOTAL  SIMULATION  TIME  =  120.0000 
CONDUCTIVITY  =  10.0000 

LATENT  HEAT  -  80,0000 

POROSITY  =  0,4000 


NODE 

NO. 

Xf  I  ) 

Yd) 

KTYPECl) 

1=SU52=SF 

3-EFFLUX 

VALUE 

ANGLE ( I ) 

1 

0.00000 

1<) ,  ooooo 

2 

0 . OOOOO 

90.00 

9 

25.00000 

10,00000 

1 

0.00000 

180.00 

3 

49.00000 

10. ooooo 

1 

o, ooooo 

1 80 . 00 

4 

49 . 90000 

10.00000 

1 

0,00000 

180.00 

5 

50.00000 

10.00000 

1 

o . ooooo 

206,57 

6 

50. 10000 

9.95000 

1 

0, ooooo 

180.00 

y 

51.00000 

9.50000 

1 

0 . ooooo 

180.00 

9 

60.00000 

5 . ooooo 

1 

0 . ooooo 

180.00 

9 

69.00000 

0,50000 

1 

0.00000 

180.00 

10 

69.90000 

0 ,  (>5000 

1 

0. ooooo 

180.00 

11 

70,00000 

o , ooooo 

1 

0.00000 

153.43 

12 

70 . ioooo 

o , ooooo 

1 

o . ooooo 

180.00 

13 

71.00000 

0.00000 

1 

O , ooooo 

180.00 

14 

95.00000 

0.00000 

1 

0 , ooooo 

180.00 

15 

120.00000 

0 ,  OOOOO 

1 

0.00000 

90.00 

16 

120.00000 

1 , ooooo 

1 

-10.00000 

90,00 

17 

95 . 00000 

l ,  ooooo 

1 

-10.00000 

180.00 

10 

71 .00000 

l , ooooo 

1 

-10.00000 

180.00 

19 

70.10000 

1 ,  OOOOO 

1 

-10.00000 

180.00 

20 

70 . OOOOO 

1,00000 

l 

-10,00000 

206.57 

21 

69.90000 

1.05000 

1 

-10, ooooo 

180.00 

22 

69.00000 

1.50000 

1 

-10.00000 

180.00 

23 

60.00000 

6,00000 

1 

- 10 . ooooo 

180.00 

24 

51 .00000 

10.50000 

1 

-10.00000 

180,00 

25 

50. 10000 

10.95000 

1 

-10.00000 

180.00 

26 

50.00000 

1 1 , ooooo 

1 

-10.00000 

153.43 

27 

49.90000 

1 1 <  ooooo 

1 

-10.00000 

180.00 

28 

49 . OOOOO 

11,00000 

1 

-10.00000 

180.00 

29 

25.00000 

l i , ooooo 

1 

-10.00000 

180.00 

30 

o , ooooo 

1 1.00000 

1 

-10.00000 

90.00 

23 


%  V  * 


New  Coordinates  of  the  Freezing  Front 
Time  =  120.0000 


Node 

X-Coord. 

Y-Coord. 

1 

0.0000 

9.6542 

2 

25 . 0000 

9.6542 

3 

49.0000 

9.6369 

4 

49.9000 

V . 5686 

5 

50.0000 

9.5203 

6 

50 . i 000 

9.5124 

7 

51 <  0000 

9.1105 

8 

60.0000 

4.6:*.S8 

9 

69.0000 

0.1376 

10 

69.9000 

-0.2377 

11 

70.0000 

-0.2365 

12 

70.1000 

-0,2807 

13 

71 , 0000 

-0.3338 

14 

95.0000 

-0,3470 

15 

120.0000 

-0.3468 

m 


* 


m 


m 


w- 

S5S 


#3 


53 


m 


„v 


»>a 
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Example  2:  Nodal  Density  and  Timestep  Size  Sensitivity  Analysis^ 

A  sensitivity  analysis  is  prepared  examining  different  time  incre¬ 
ments  and  nodal  point  densities  and  the  resulting  effects  on  CVBFR1 
modeling  results.  Figure  3.1  shows  the  different  nodal  densities  and 
Fig.  3.2  shows  the  results  from  the  several  CVBEM  models.  From  the  analysis 
it  appears  that  a  small  timestep  (6-hours)  is  preferred,  but  a  large 
timestep  such  as  60  hours  results  in  a  relative  error  with  respect  to 
the  one-dimensional  Stefan  solution  of  only  2  percent.  Additionally,  a 
relatively  sparse  nodal  density  of  only  30  nodes  results  in  a  satis¬ 
factory  approximation. 

Example  3:  Comparison  to  Two-Dimensional  Domain  Model inq  Results 


The  CVBFR1  modeling  results  for  the  previous  example  are  compared  to 


results  from  a  Nodal  Domain  Integration  (NDI)  two-dimensional  phase  change 
model  in  Fig.  3.3.  The  NDI  model  is  based  upon  an  isothermal  soil -water 
phase  change  approximation,  and  uses  an  apparent  heat  capacity  approach 
to  model  the  freezing  front  evolution  in  the  fixed  grid  domain  model. 


»j 

it 


'i 
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Fig.  3.3b.  Comparison  of  CV6EM  and  NDI  Modeling  Results 
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APPENDIX  A: 

COMPLEX  VARIABLE  BOUNDARY  ELEMENT  METHOD 


Hromadka  and  Guymon  (1984c)  present  a  detailed  development  of  the  CVBEM. 

A  comprehensive  presentation  of  the  method  is  given  in  Hromadka,  (1984,  1987). 

A  feature  available  with  the  CVBEM  is  the  generation  of  a  relative  error  meas¬ 
ure  which  can  be  used  to  match  the  known  boundary  condition  values  of  the  prob¬ 
lem.  Consequently,  the  method  can  be  used  to  develop  a  highly  accurate  approx¬ 
imation  function  for  the  Laplace  equation  and  yet  provide  a  descriptive  rela¬ 
tive  error  distribution  for  analysis  purposes.  Because  the  main  objective  of 
this  paper  is  to  analyze  the  numerical  error  in  solving  (5),  it  is  noted  that 
the  Laplace  equation  is  solved  throughout  the  problem  domain  (if  homogeneous) 
or  in  connected  subregions  (if  inhomogeneous).  Many  anisotropic  effects  can  be 
accomodated  by  the  usual  rescaling  procedures  or  by  subdividing  the  total  do¬ 
main  into  easier-to-handle  subproblems.  The  CVBEM  is  then  applied  to  the  prob¬ 
lem  domain(s)  as  discussed  in  the  following. 

Let  Q  be  a  simply  connected  domain  with  boundary  r  where  r  Is  a  simple 
closed  contour  (Fig.  A1 ) .  Discretize  r  by  m  nodal  points  into  m  boundary 
elements  such  that  a  node  Is  placed  at  every  angle  point  on  r  (Fig.  A2).  Each 
boundary  element  is  defined  by 

Tj  »  {z:  z  »z(s)  where  z(s)  ■  Zj  +  (Zj+1  -  Zj)s,  0  <s  s  1),  jfm  (Al) 
with  the  exception  that  on  the  last  element, 

r  *  (z:  z  * z( s )  where  z(s)  *  zm  +(z.  -z Js,  0  <  s  <  1  > 

Hi  hi  1  m 


r'  M  rj 

j-i 


C-  C'aV 


•  ^  O'.  O  v'  »  *  - 


Let  each  Tj  be  discretized  by  (k+1)  evenly  spaced  nodes  (k  *1)  such  that 
Is  subdivided  Into  k  equi length  segments  (Fig.  A3).  Then  Tj  Is  said  to  be 
a  (k+l)-node  element.  Prom  Fig.  A3,  each  Tj  has  an  associated  nodal  coordinate 
system  such  that  zja  -  and  zjfk+1  ■  tj.j  -  zJtl>r 
On  each  r^,  define  a  local  coordinate  system  by 


*  ♦  <V1  -  *j)s 


(A3) 


.here  dCj  ■  (zJ  k+1  -  ZjJds. 

On  each  (k+l)-node  element  1^,  a  set  of  order  k  polynomial  basis  functions 
are  uniquely  defined  by 


(s)  •  a 


J.1.0 


♦  a 


J.1,1 


s 


+**,+aJ,1.k  s 


(A4) 


where  1  «  1,2, •••  ,(k+l)  and  0  <s  <  1,  and  where 


j  »k 


zj.n  *zj,l 

s  , 

zj,k+l  ”zj,l 

1,  n  *  1 
O.nM 


(A5) 


The  basis  functions  are  further  defined  to  have  the  property  that  for  z,  er 


j  .1 


v  ZJ  ,k+l  " zj  ,1 


;  -  z 


j.i 


Zj ,k+l  ’zj,l  > 


.  ur, 


(A6) 


o  .  <tr. 


Let  w(z)  be  analytic  on  Q  U r.  That  is,  let  w(z)  be  the  solution  (unknown) 
to  the  steady-state  boundary  condition  problem  being  considered.  At 
each  nodal  point  on  r,  define  a  specified  nodal  value  by  (Fig.  A3) 


LEGEND 


•  ELEMENT  ENDNODE 
o  ELEMENT  INTERIOR  NODE 

FIG.  A3.  ( k*  I)-  NODE  BOUNDARY  ELEMENT  Ij  NODAL  DEFINITIONS 


FIG.  A4.  BRANCH-CUT  OF  LN  (z-  f)  FUNCTION,  feT 


where  from  Fig.  A3,  “j  \  *  wj  *  “j-l.MT 


Using  (A6)  and  (A7),  an  order  k  global  trial  function  is  defined  by 


(A8) 


From  ( A8) ,  the  global  trial  function  is  continuous  on  r.  An  approximation 
function  <2^(z)  (Hromadka,  1984,  1987)  is  defined  by  the  Cauchy  integral 


^(z)  -  — 

K  '  2ir1 


f  GkU)  dc 


,  Zell,  z{  r 


(A9) 


c  -2 


Because  the  derivative  of  ^(z)  exists  for  all  z  eft,  then  u^U)  Is  analytic 
In  ft  and  exactly  solves  the  Laplace  equation  In  ft. 

Expanding  (A9)  and  using  (A2)  gives 


’  Gk(?)  dc 
C  -  2 


m 


j-1 


GkU)  dc 
c  -  2 


Integrating  on  boundary  element  j  gives  (Hromadka,  1984,  1987) 


(A10) 


'  Gk(c)  dc 

c  -  1 


*  Rj"’ul  * 
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z  -  z 


(All) 


k- 1 

where  Rj  (z)  is  an  order  (k-1)  complex  polynomial  resulting  from  the  circuit 

around  point  z  (see  Fig.  A4)  and  y,  is  equal  to  (z-z.)/(z.  .  -  z  .)• 

J  J  J 

Thus,  the  CVBEM  results  in  the  approximation  function 


^  M  l  KW(1)  *  ?  “J.t  "j>j> 


z  -  z 


J+l 


(A12) 


or  In  a  simpler  form  (Hromadka,  1984,  1987) 

<1  (z)  -  R k(z)  ♦  —  I  in  (2 -zA  l  T  k  (A13) 

^  2ir1  j  J  1  1 

where  *  iL^  ^  ^  <Yj_ j )  -  Wj  ^  (yj).  and  R*(z)  follows  from 

(A12). 

The  approximation  function  of  (A13)  exactly  satisfies  the  governing  flow 
equation  in  the  problem  domain  (1  for  the  approximated  boundary  conditions  on 
the  problem  boundary,  r.  Because  ^(z)  Is  analytic  on  Q,  then  the  maximum 
relative  error  of  |u(z)  -  ^ (z ) |  must  occur  on  r.  Consequently,  the  total 
approximation  error  can  be  simply  evaluated  on  r  with  the  corresponding  errors 
in  the  interior  of  ft  being  less  in  magnitude.  Because  the  boundary  conditions 
used  to  evaluate  (A1 3)  are  known  continuously  on  r,  then  uL(z)  can  be  determined 
within  arbitrary  accuracy  by  the  addition  of  nodal  points  on  r  due  to  (without 


proof) 


lim 


G*U)  d< 


2ni  lim  <1  ( z )  = 


maxlr,  |-*o 
_ _ _ 


f  wU)  dc 


=  2tt1  uj(z)  (A14) 


•  v . 


APPENDIX  B: 


THE  APPROXIMATIVE  BOUNDARY 
FOR  CVBEM  ERROR  ANALYSIS 

Generally,  the  prescribed  boundary  conditions  are  values  of 
constant  $  or  0  on  each  Ij.  These  values  correspond  to  level  curves 
of  the  analytic  function  ui(z).  *  0  ♦  i0.  After  determining  a  u(z),  it 
is  convenient  to  determine  an  approximative  boundary  f  which  corresponds 
to  the  level  curves  of  <j(.z)  »  0  ♦  10  which  are  specified  as  the  prescribed 

boundary  conditions.  The  resulting  contour  f  is  a  visual 
representation  of  approximation  error,  and  r  coincident  with  V  implies 
that  u(z)  *  uj(z ) .  Additional  collocation  points  are  located  at  regions 
where  ?  deviates  substantially  from  r. 

A  difficulty  in  using  this  method  of  locating  collocation  points 
Is  that  the  contour  r  cannot  be  determined  for  points  z  outside  of  flUT. 

To  proceed,  an  analytic  continuation  of  d>( z )  to  the  exterior  is  achieved 
by  rewriting  the  integral  function  (A9)  in  terms  of 

1  GU)d<;  m 

— -  -  »  R  (z)  +  l  (a,  +  1 S  . )  ( Z  -  Z  . )  Ln  (z-zj  (Bl) 

2-rri  J  c-z  J-l  J  J  J  J 

* 

where  aj  and  8j  are  real  numbers;  and  Ln  (z-Zj)  is  a  principle  value 
logarithm  with  branch-cuts  drawn  normal  to  r  from  each  branch  point  Zj 
such  as  shown  in  Fig.  Bl .  The  resulting  approximation  is  analytic  every¬ 
where  except  on  each  branch-cut.  The  Rt(z)  function  In  Eq.  (Bl)  Is  a 
first  order  reference  polynomial  which  results  due  to  the  Integration 
circuit  of  2 it  radians  along  r.  If  u(z)  Is  not  a  first  order  polynomial. 
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Implementation  on  a  computer  Is  direct  although  considerable  com* 
pu tat Ion  effort  Is  required.  One  strategy  for  using  this  technique  Is 
to  subdivide  each  Tj  with  several  Internal  points  (about  4  to  6)  and 
determine  u(z)  at  each  point.  Next,  ,r  Is  located  by  a  Newton- Raphson 
stepping  procedure  In  locating  where  w(z)  matches  the  prescribed  level 
curve.  Thus,  several  evaluations  of  w(z)  are  needed  to  locate  a  single 
?  point.  The  end  product,  however,  may  be  considered  very  useful  since 
It  can  be  argued  that  u>(z)  Is  the  exact  solution  to  the  boundary  value 

<<N  A 

problem  with  r  transformed  to  r,  and  r  Is  a  visual  Indication  of 
approximation  error. 

The  use  of  the  method  discussed  for  locating  additional  collocation 
points  on  r  Is  demonstrated  by  application  of  the  CVBEM  for  solving  2  steady 

"  V 

state  heat  transfer  problems.  The  problems  considered  each  Involve  a  different^ 
geometry  and  set  of  boundary  conditions  of  the  Dlrichlet  class.  The  analytic 
solution  to  the  problems  are  Included  In  F1g.B2,  Each  solution  satisfies 

the  Laplace  equation  and  Is  defined  as  a  function  of  a  local  coordinate  x-y  < 

system  with  an  origin  specified  as  shown  In  the  figures.  On  the  problem 

boundaries,  r,  the  potential  function  or  temperature  Is  also  a  continuous  T; 

function  of  position  defined  by 


♦(»  c  D  -  £  (*2  ♦  y2) 


From  (B2),  It  Is  seen  that  the  boundary  conditions  are  not  level  curves; 
consequently,  the  determination  of  an  approximative  boundary  f 
requires  further  definition.  In  these  applications,  the  problem  Is  approached %•'.%- 
by  using  the  statement 


T  i  •  z  :  *U)  ■  |  (x*  ♦  yJ)  *  £  |  z  |  • 

2  4 


.v-vvv.  .* 


.v.vj /.v/; 


The  strategy  of  working  with  level  curves  (l.e.  <p  •  for  1  e  ij,  J  »  1,2, •• 
follows  analogously. 

The  two  applications  Illustrate  the  development  of  CVBEI1  approximation 
functions  which  exactly  satisfy  the  governing  partial  differential  equation 
(Laplace  equation)  in  fl  and  approximately  satisfy  the  boundary  conditions 
which  are  continuously  specified  on  r.  The  subsequent  figures  illustrate  the 
CVBEM  error  evaluations  along  r  for  evenly  spaced  nodal  placements 
for  each  problem  boundary. 
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B4  APPROXIMATIVE  BOUNDARIES  FOR  FIVE  NODAL  POINT 

DISTRIBUTIONS 


n  n  o »  o  oo  -j  o  o  o  n  n  r>  o  o  .->  n  n  o  n  n  n  o  n  o 
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APPENDIX  C:  PROGRAM  CVBFR1  Fortran  Listing 


MAIN  PROGRAM 

THIS  CAUCHY  PROGRAM  (  FREEZING  OR  THAWING  FRONT  ADVANCEMENT  > 

USES  SUBROUTINES  CAUCH1 . CAUC H2 . CAUCH3 . CAUCH4 , CAUCH5 . HQM . ANG . FRONT 

BASED  ON  THE  APPROXIMATION  FUNCTION 

IMPLICIT  DOUBLE  PRECISION (A-H.O-Z) 

COMMQN/BLK  1/X(100> 

COMMON/BLK  2/Y<100> 

COMMON/BLK  3/KT YPE ( 100 > 

COMMON/BLK  4/VALUE  <  1 00  ) 

VIRTUAL  P(IOO.IOO) 

COMMON/BLK  6/S<100) 

COMMON/BLK  7/ANGLE (100' 

COMMON/BLK  9/NAT ( 100 ) 

DIMENSION  REX< 100) . REY< 100) 

DIMENSION  HIY(IOO) 

OPEN  DATA  FILES 

WRD*1 
NWT  =  0 

OPEN(UNIT=NRD.NAME= 'CAUFRT , DAT ' ? TYPE* ' OLD ' > 
OPEN<UNIT=NUT»NAME='CAUCHY.ANS' . TYPE='NEW' > 

READ  DATA 

...note:  node  number  plus  number  of  efflux  bc  , 

(NNODP=NNOD+NNAT)  CAN  NOT  EXCEED  *  100* 

READ ( NRD  » * ) KODE 
READ (NRD. » ) NNOD . NFRS . NFRE 
READ ( NRD .  * > COND » XL AT . POR 
READ (NRD. * > BEL T . S IMUL . OUT . ID 
WR I TE ( NWT -601) BELT  • S I MUL . COND . XL AT . POR 
01  FORMAT <///.6X.  'TIME  INCREMENT  =  '  .  F8  =  4 .  /  >  6X  .  '  TOT  AL  SIMULATION', 
i'  time  =  ' .F8.4./.6X. 'CONDUCTIUITY  =  '• F3  =  4 •  /, 6X LATENT  '  . 

2'  HEAT  =  ' .F8. 4./. 6X. 'POROSITY  =  '.F6.4./> 

.  ,  . "ALUE  OF  EFFLUX  B.C  =  EFFLUX/CONDUCTIUITY 
DO  7  1  =  1. NNOD 

READ (NRD. * ) 'X< I ? . Y< I > >KTYPE( I > . VALVE ' 1 ) > 

CALL  ANG ( NNOD  > 

WRITER  NWT  » 1 0  > 

0  FORMAT ( 6X.  'NODE' » 6X . ' X ( I ) ' . 6X ? ' Y ( I > ' . 4X .  'KTYPE (  I  >  '  .  3X .  'VALUE  '  . 
15X. ' ANGLE ( I ) ' ./.  7X. 'NO .  ' . 24X. ' 1=SV. 2  =  SF' ./ • 35X. ' 3=EFFLUX  '  > 

DO  9  1*1. NNOD 

WRITE (NUT. 8) I »X< I ) . Y ( I ) . KTYPE <  I  > .VALUE' I > . ANGLE  (  I  > 

FORMAT (3X. 15 . 5X . 2F 10 . 5 . 15 . 5X . F 10 = 5 . F 10 . 2 ) 

CONTINUE 

WR I TE ( NWT . 602  > 

02  FORMAT ( 72 (  ' - '  >  ) 

CHECK  NATURAL  OR  EFFLUX  BOUNDARY  CONDITION 
WNAT =  0 

DO  3  1=1. NNOD 
X  < I >  =X  <  I  > 

Y  (  I )  =  Y  (  I  ) 

IF  <KTYPE( I > i NE . 3 > GO  TO  3 
NNAT  =  NNAT  6 1 


BK 


non 


NNQDP=NNOD+NNAT 
NAT ( I ) =NNQDP 
CONTINUE 

IF'NNAT, EQ.0>NN0DP=NN0D 

PREPARE  GLOBAL  MATRICES 

,  . ZERO  ARRAYS 

ITER*IFIX<  SIMUL/DELT ) 

IOUT*IFIX(OUT/DELT> 

KOUT  =0 

DO  9999  1 1 1 1*1 ? ITER 
KOUT  =  KOUT  +  1 
DO  5  1=1* NNODP 
S(D*0. 

DO  6  1  =  1*  NNODP 
DO  6  11*1*  NNODP 
P(  I  *  II  )=0  , 

DO  1000  J= 1  *  NNOD 
..ACCOMODATE  DIAGONAL  NODE 
I  =  J  - 1 

IF  <  I  .  EQ  •  0 ) I =NNOD 
K=  J  +  l 

IF  (  K  :  GT  :  NNOD  K  =  1 

CALL  CAUCHK  J  *  I  ?  K  *  A  *  B  *  C  *  D  > 

A  J  =  A 

BJ=ANGLE(J> /ISO. *3.141593 

CALL  CAUCH2<  J*I*K* A*B*C* D* AJfBJ* P' 

,  , ACCOMODATE  REMAINING  CONTOUR  NODAL  POINTS 
NELE=NNGD-2 
DO  500  K= 1 • NELE 
M  =  J  +  K 

IF  (  M  <  GT  NNOD ) M  =  M-NNOD 
N  =  M  + 1 

IF (N .GT : NNOD ' N=N-NNOP 
CALL  CAUCHK  J*M*N*  A*B*C*D> 

CALL  CAUCH2C J * M * N *  A  * B * C * D *  A J * B J * P > 

0  CONTINUE 
00  CONTINUE 

PREPARE  RELATIVE  ERROR  ANALYSIS 

CALL  CAUCH3 ( NNODP  *  NWT  *  P  > 

TIME  =  DELT*FLOAT  < II II > 

IF'KOUT  EQ.  I  OUT > CALL  CAUCH4 ( NNOD * NWT • T I  ME ? I D > 

ASSIGN  BOUNDARY  NODAL  POINT  VALUES 

DO  7010  1  =  1  * NNOD 
IF(KTYPE<I) . EQ . 2  >  GO  TO  7015 
IF(KTYPE( I > = EQ:3>G0  TO  7016 
REX(I>=VALUE'I> 

RE Y  < I >  *5  <  I  > 

GOTO  7010 
15  REX ( I >  =S < I ) 

RE Y ( I ) = VALUE ( I  > 

GOTO  ^010 
'16  1 1  *NAT  (  I ) 

REX ( I ) *S  < I ) 

RE Y ( I ) =S ( II  > 

HO  CONTINUE 


5 


kV 


-  •F  -V-V -  *£**./£ 


CALCULATE  K  £  L  A  T  I  u  r  frrqr  “ALUES 


CALL  HOM'REX.REY.NNOD.NWT.KOUT. I0UT.H1Y.ID) 

UPDATE  THE  NEW  INTERNAL  ANGLES  AND  POSITIONS 
OF  THE  FREEZING  OR  THAWING  FRONT 

CALL  FRONT < NNOD » NFRS » NFRE » COND « XL AT . POR . DEL T . HI Y.KQDE) 

OUTPUT  THE  NEW  POSITIONS  OF  THE  FREEZING  OR  THAWING  FRONT 

IF'KOUT  NE:  1 0  U  T  >  G  0  T0  999? 

IF'ID  .NE:  0>WRITE<NWT. 6031TIME 
FORMAT </*6X. 'TIME  =  '.F9.2./> 

WR I TE ( NMT  » 603  > 

FORMAT < /• 6X . 'NEW  COORDINATES  OF  THE  FREEZING  OP  THAWING  FRONT 
1 . .  12X.  '  NODE  '  .  4X  t  '  X -COORD  :  '  .3X.  '  Y -COORD  -.'•/> 

DO  400  I =NFRS  > NFRE 
X  A  =  X  <  I  > 

Y  A  =  Y  (  I  > 

WRITE'NUT.604>I.XA. YA 
FORMAT' 1 IX. I3.5X.F9 : 4.6X.F8 : 4) 

CONTINUE 
WRITER  NWT . 402  > 

K  0  U  T  =  0 
°  CONTINUE 

CLOSE(UNITsNRD) 

CLOSE(UNIT =NWT > 


Cl  Cl  u 


SUBROUTINE  CAUCHi 


SUBROUTINE  CAUCHI  (JfM-N-AfB.-CrP' 

C  IMPLICIT  DOUBLE  PPEC I S I  ON < A -H » 0- 2 > 

COMMON/BLK  1/XUOO' 

COMMON/BLK  2/Y<100) 

COMMON/BLK  3/KTYFE'100> 

COMMON/BLK  7/ANGLE < 100 ) 

C 

C  SUBROUTINE  T0  DETERMINE  BOUNDARY  ELEMENT  GEOMETRIC  VALUES 
r 

C ..  .CALCULATE  VECTOR  LENGTHS 

XLN*SQRT< < Y(N)-Y( J) > **2+ ' X < N > -X <  J) >**2> 

XLM*SGRT  < (Y<M)-Y< J> > **2+ < X < M > -X < J) >**2> 

XXX=XLN/XLM 
A  =  ALOG ( XXX  > 

C  A=DLOG( XLN/XLM) 

C .. .DETERMINE  ANGLE  ARITHMETIC 
ZMX  = ( X ( M  > -X ( J ) ) / X L M 
ZMY*< Y<M)-Y< J) ) / X  L  M 
ZNX = ( X ( N ) -X ( J ) )/XLN 
ZN Y=  <  Y ( N  '  - Y ( J> > /XLN 
CALL  CAUCH5I ZNX , ZNY » ANGLEN) 

CALL  CAUCH5<ZMX»ZMY» ANGLEM) 

B=ANGLEN-ANGLEM 

C. . .ACCOMODATE  CENTRAL  ANGLE  DETERMINATION  BEING  BACKWARDS 
I F  <  M i E  Q :  <  J  —  1 >  .OR:  N : EQ . <  J+l ) ) GO  TO  98 

c..  accomodate  branch-cut  effects 

IF ( B . LT ; -3. 14159)8=8+6.2831853 
IF  <  B . GT .  3. 14159>B  =  B-6. 2831853 

GOTO  90 

08  CONTINUE 

D  =  A  N  G  L  E ( J) 
oo  CONTINUE 

r 

C.. COMPLEX  VARIABLE  ARITHMETIC 

C 

F*OC<N>-X'M)>**2+'Y'N>-Y<M>>*»2 

C=A»(X(N)-X(M))-B*(Y(M)-Y(N)> 

D=B*<X'N'-X<M))+A*<Y'M)-Y<N'> 

C  =  C  /  F 
D  =  D/F 
RETURN 
END 


o  it  n  o  it  n  nn  n  n  o  o  o  o  non  o 


SUBROUTINE  CAUCH2 


SUBROUTINE  CAUCH2 < J r M » N» A » B r C f D r A J» BJ » P > 

IMPLICIT  DOUBLE  PREC I S I  ON ( A-H » 0-Z ) 

COMMON/BLK  1/X(100) 

COMMON/BLK  2/Y(100> 

COMMON/BLK  3/KTYPE<100> 

COMMON/BLK  4/VALUE<100> 

VIRTUAL  P(100fl00> 

COMMON/BLK  6/S' 100) 

COMMON/BLK  7/ANGLE(lOO> 

COMMON/BLK  8/NATdOO) 

SUBROUINE  TO  ASSEMBLE  BOUNDARY  ELEMENTS 

INTO  GLOBAL  MATRIX  *P*  WITH  VECTOR  OF  CONST  ANTS  r ‘S' 


SUBROUINE  TO  ASSEMBLE  BOUNDARY  ELEMENTS 

INTO  GLOBAL  MATRIX  *P*  WITH  VECTOR  OF  CONS TANTS • * S * 

F=AJ*AJ+BJ*BJ 
AZ=-A J/F 
BZ=-BJ/F 
JJ= J-l 

IF  <  M . EQ . J J  >  GOTO  100 
JJ=JF1 

IF'N  EQ>JJ)GOTQ  100 
...ELEMENT  DOES  NOT  CONTAIN  NODE  *J* 

C1  =  'X'J)-X'M>>*C-'.  Y'J)-Y'M>)*B  +  1. 

C  2  =  '  X 'J)-X'M>'*D+'Y(J)-Y'M>>*C 
C3*'X<J>-X'N>>*C-'Y<J)-Y<N>>*B+1. 

C4= ( X '  J'-X'N)  >*D+'Y(  J  >  -  Y  (  N  >  >tC 

CCl=Clf AZ-BZ4C2 

CC2=C1*BZFC2*AZ 

CC3=C3*AZ-C4fBZ 

0C4=C4*AZ*BZ<C3 

C1=CC1 

r  ">  =  c  r  t 

C3=CC3 

C  4  *  C  C  4 

...ASSIGN  COEFFICIENTS  to  UNKNOWN  HARMONIC  VARIABLE 
I F ' K  TYPE ' J  > . EQ : 1 ) GO  TO  5 

..DIAGONAL  NODAL  UNKNOWN  HARMONIC  IS  THE  STATE  VARIABLE 
. . . USE  REAL  EQUATION 
G 1 =-C  3 
G2*C4 
G3*C  1 
G4*-C2 
GO  TO  9 

..DIAGONAL  UNKNOWN  HARMONIC  IS  THE  STREAM  FUNCTION 
...USE  IMAGINARY  EQUATION 
0 1 =-C  4 
G2=-C3 
G  3  =  C  2 
G4-C1 

9  IF<KTYPE<M)  EQ.2)G0T0  10 

IF  <  KTYPE  '  M  )  :  EQ  =  3)G0T0  15 
C... STATE  VARIABLE  SPECIFIED  FOR  NODE  'M* 

S' J'=S<  J)-(G1 >tVALUE<M> 

P' J»M'=P' J.M)+'G2) 


C>  C.I  Cl  Cl  Cl 


X 


..EFFLUX  SPECIFIED  ^OR  ^ODE  *  M  * 

S( J)*S( J> 

P< J»H>»P'  J  •  M  '  +  G 1 
MF  =  NA  T (  M  ' 

P( J  *  MF ) =P  <  J  »  MF ) +G2 
QO  TO  50 

..STREAM  FUNCTION  SPECIFIED  FOR  NODE  *m* 

S< J)=S< J)-(G2)*VALUE<M> 

P< J»M) =  P< J»M)+(G1 ) 

IF  < KTYPE ( N ) . EQ . 2 ) GOTO  60 
IF'KTYPE'.N)  EQ.3>G0T0  45 
..STATE  VARIABLE  SPECIFIED  FOR  NODE  *N* 

S'  J>=S< J>-'G3>1tVALUE'N> 

P< JrN>»P< JfN)+(G4> 

GO  TO  250 

..EFFLUX  SPECIFIED  FOR  NODE  *N* 

S  (  J  >  =S  <  J  > 

P<  J»N>*P'  J.*N)+G3 
N  F  =  N  A  T  <  N  > 

P( J»NF>=P( J » NF  >  +G4 
GO  TO  250 

..STREAM  FUNCTION  SPECIFIED  FOR  NODE  *N* 

S  J  >  =S  <  J  >  -  '  G4  >  ♦ VALUE  (  N ) 
p  <  j  .  N )  =  P  (  j .  N  >  -f  <  G  3  > 

GO  TO  250 

BOUNDARY  ELEMENT  CONTAINS  NODE  ‘J1 

l  IF ' NTYPE ' J > . EQ  .  2  .OR.  KT YPE ' J )  EO  3 > GO  TO  110 
.STATE  VARIABLE  SPECIFIED  FOR  NODE  * J* 

USE  IMAGINARY  EQUATION 

IF  (KTYPE(N)  .  EQ  .  1  )  P  (  J  f  N  )  =P  (  J-N.'-FAZ 

IF  < KTYPE  ' N I  . EQ  .  1 >  S<  J ) =S<  J ? -BZTVALUE ( N  > 

IF  'KTYPE'  N  ,  EQ.  2>P<  Jf  N)=P(  JrN'+BZ 
IF  <  KTYPE  <  N  '  .  EQ  .  2  >  S  <  J  )  *S  (  J  >  -AZf VALUE  '  N  > 
IF'KTYPE'N) . NE . 3>G0  TO  113 
.EFFLUX  SPECIFIED  FOR  NODE  *N* 

S  (  J  '  =  S  '  J  ) 

P  '  J  •  N  =  P  <  J  •  N  )  +  B  Z 
N  F  =  N  A  T  '  N  > 

p  <  J  •  N  F  '  =  R  i  j  .  n  F  >  *  A  Z 
IF  ' K  TYPE ' M  >  .EQ.  2  .'GOTO  115 
IF  i.  KTYPE  <  M  >  :  EQ  .  3  >GOTO  114 
S' J)=S( J)fBZ*VALUE(M> 

P<  J»M)*P( J?M)-AZ 
GO  TO  200 

S( J)=S' J)+AZ*VALUE'M> 

P< JfM)=P( J»M)-BZ 
GO  TO  200 

: EFFLUX  SPECIFIED  FOR  NODE  *  M ■ 

S  '  J  >  =  S  <  J  > 

P( J-M)=P( J  »  M ) -BZ 
M  F  =  N  A  T  '  M  > 

c'(  J-MF)=P(  J  ?  MF  >  -  AZ 
GO  TO  200 

.STREAM  FUNCTION  SPECIFIED  FOR  NODE  ’J1 
IF (KTYPE(N) .NE. 1)G0T0  120 
S<  J>=S<  J>-AZ*VALUE(N> 

P< JrN)=P( Jf N)-BZ 
GO  TO  130 

IF ' KTYPE  <  N  > .NE . 3  >  GO  TO  111 


:  .  .  ,  EF-LUX  3 F £ C  I  F  I  E I'  FOP  NODE  1  n  • 

S  (  J  ?  =  S  (  J  > 

P(  J,N>=P<  J»N>+AZ 
NF  =  NAT < N ) 

P ( J  »  NF  >  =P  <  J  f  NF  > -PZ 
GO  TO  130 

111  S<  J)=S(  J)+BZ#VAL'JE'N' 

P<  J»N>*P<  J*N)+AZ 

130  IF(KTYPE(H) .NE : 1 >QQ  TO  140 
S<  J)=*S<  J ) +AZ<VALUE ( H  > 

P< J,M>=P( J>M)PBZ 
GO  TO  200 

140  IF  (  KTYPE  (  M  )  •.  NE-;  3  ?  GO  TO  112 

C... EFFLUX  SPECIFIED  FOR  NODE  ‘N* 
S(J>*S<J> 

P< Jf H)=P( J»M)-AZ 
NF  =  NAT  < M > 

F< JrHF)=P( Jr HFl+PZ 
GO  T0  200 

112  SC  J  '  =  S  C  J>-BZtUALUECM> 

P( J,H>=PC J  •  M  >  -  A  Z 

200  IF  C  KTYPE  C J) .  EQ  : 3) GO  TO  100 
P' J. J)=p< J. J  > - 1  , 

GO  TO  200 

r,,. EFFLUX  SPECIFICED  FOP  NODE  *J* 

100  JF  =  NAT  < J  > 

*F=NATCM> 

PZZ=CX(J>-X(N' >  4  *  2  +  C  Y  C J  >  —  Y  <  m  >  1**2 
DZ2  =  SQPT(,PZZ' 
c  /  jp  \  =s  (  jF  '  -I'ALUE  <  J  i  KDZZ 
P  <  JF  •  JF  >  =  1  . 

IF ( KTYPE  <  M  > : NE  3 ' P <  JF  •  N > = -  1  , 

IF i KTYPE <  N > : EQ . 3>P( JF  t MFl =-  1  . 

P  <  j  .  j  >  =  p  t  j  .  J  )  -  1  , 

200  CONTINUE 
RETURN 
END 
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C  SUBROUTINE  CAUCH3 

- - 

SUBROUTINE  CAUCH3  <  NNOD • NUT • P ) 

IMPLICIT  DOUBLE  PREC I S I  ON < A-H * O-Z  ' 

THIS  SUBROUTINE  SOLVES  A  NNOOxNNOD  MATRIX  SYSTEM 
GAUSSIAN  ELIMINATION  METHOD  USED. 

VIRTUAL  P  ( 100? 100) 

COMMON /BLK  4/S<100> 

N 1 =NNOD- 1 
DO  100  K=1*N1 
Kl=K-tl 
C  =  P  (  K  f  K  > 

IF  < APS  <  C  )  -  000001 >10? 10*70 
DO  20  J  =  K  i r  N N 0 D 

I F ( A B S ( P ( J  •  K ) > - . 000001  >  2  0  f  2  0  ? 15 
DO  14  L=K • NNOD 
r  =  p  r  k;  •  L  > 

P  (  K  •  L  >  =P  (  J  .•  L  > 

P  I  j ,  l  >  =  C 

C  =  s  f N ' 

5  (  K  >  =  S  (  J  > 
g  <  j  i  =r 
r  =  p  <  k;  .  k;  > 

GO  Tn  ->n 
CONTINUE 

WRITE' NUT.  i  >k; 

FORMAT'  1X>  ' "  "SINGULARITY  IN  ROW  '  •  I  5  > 
r.n  t  g  300 
r  =  p  t K . k 1 
DO  90  J=K 1 • NNOD 
p  t  (s; .  j  )  =  P  t  x  ,  j  >  /  C 

S  (  K  )  “  S  <.  K. )  /  C 
DO  °0  I  =M  •  NNOD 
C  =  P< I f  K  > 
n  n  00  J  =  K 1  •  N  N  0  D 
P < I • J  ) =P < I • J  > -C#P ( K • J > 

S  I  )  =  S  (  I  >  -  C  f  S  (  K  ) 

CONTINUE 

I F  '  ABS  P  NNOD  •  NNOD  >  '  -  000001 '30-30*  120 
S  1.  NNOD  )  =S  <  NNOD) /P  1  NNOD  •  NNOD  > 

DO  200  L  =  1 •  N  1 
K=NNOD-L 
K  1  =  K  + 1 

DO  200  J  =  K  1 >  N  N  0  D 
S  ( K  )  *S  ( K. )  -P  <  K  f  J  >  tS  J  > 

CONTINUE 
RETURN 
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SUBROUTINE  CAUCH4 


SUBROUTINE  CAUCH4  '  N.N0D  .•  NWT  ?  TIME  .•  ID  > 
IMPLICIT  DOUBLE  PREC I S ION ' A-H • 0-2 > 
CQMMQN/BLK  3/KTYPE'.100> 

COMMON/BLK  4/VALUE  <100) 

COMMON/BLK  6/S(100) 

COMMON/BLK  8/NAT' 100' 


SUBROUTINE  FOR  OUTPUT 

IF (IB  N£ ;  0)  RETURN 
URITE(NWTflO)  TIME 

FORMAT  (////// f  40.X  >  'CAUCHY  PROGRAM  RESULTS  '•/•  4X  T  IME  = 
WRITE (NWT? 12 > 

FORMAT  (/i-AX?  'NODE'  •  6X  •  '  ST  ATE  '  1 1 4X »  '  STREAM  '  •/•  5X  •  '  NUMBER  ' 
C3X  • 'VARIABLE' ■ 12X? 'FUNCTION' > 

DO  SO  1  =  1  •  N N 0 D 
IF(KTYPE(I)  .NE.DOO  TO  20 
1 1  =  N  A  T  (  I  ) 

WR  I  TE  (  NWT  .•  55  >  I  .•  S  (  I  >  •  S  (  1 1  > 

IF  (KTYPE(  I  )  .  EQ  :  1  )URITE(NWTf  53)  I  f  VALUE'.  I  >  ?  S(  I  > 

IF  (  KT YPE  < I > : EG : 2 > WRITE ' NWT ? 35 > I » S < I >  ? VALUE'  I ) 

FORMAT!  3Xf  15  fSXfFlO.  4  »  10X  » FI 0  .•  4  > 

CONTINUE 

RETURN 

END 
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C  SUBROUTINE  CAUCH5 

C - 

SUBROUTINE  C AUCH5 ' X , Y f ANGLE ) 
IMPLICIT  DOUBLE  PREC I S I  ON < A-H • 0- Z > 


THIS  SUBROUTINE  DETERMINES  THE  POS’TIUE  ANGLE 
C  OF  COMPLEX  POINT  X+iY  WITH  RESPECT  TO  THE  ORIGIN 

C 

P I  =  ACOS  < - 1  :  > 

I F  <  X  .  E  Q  .  0  .  : AND  >  Y . GT , 0 = > ANGLE* . 5*P I 

IF'X.EQ.O.  .AND.  Y.LT.O. >ANGLE=1.5*PI 
I F  <  X . GT  .  0  :  .AND.  Y  .  GE  .  0  ,  )  ANGLE  =  ATAN  <  Y/X  ) 

I F ( X . LT . 0  j  .AND.  Y.GE.O. )ANGLE  =  PI-ATAN(-Y/X> 

I F  <  X . L  T , 0 .  .AND,  Y  .  LT  ,  0 .  > ANGLE=P I  +  ATAN ( Y/X > 

I F  ■:  X  .  G  T  0  .AND.  Y  .  L  T  .  0  .  >  ANGLE  =  2  .  *P  I  -  AT  AN  '  -  Y /X  > 
RETURN 


SUBROUTINE  NOM 


SUBROUTINE  HOM  <  REX  f  REY  ?  NNOD  ?  NUT  f  KOL'T  •  I  OUT  •  H 1 Y  •  I P  ) 

IMPLICIT  DOUBLE  PRECISION ( A-H » O-Z ) 

THIS  SUBROUTINE  CALCULATES  THE  LIMITING  NODAL  POINT  VALUES 
OF  THE  ANALYTIC  HI  APPROXIMATION  FUNCTION 

CCMMON/BLK  1/X'lOO) 

CQMMON/BLK  2/Y(100) 

COMMON/ BLK  ?/ANGLE  <  1 00  ) 

DIMENSION  H1X< 100) >H1Y< 100) 

DIMENSION  REX  < 100  > » REY (100) 

MAIN  LOOP 


KPsKOUT 

IF'.IB  NE .  O)K0UT=9999 
DO  20  J  =  1  •  N N 0 D 
H1X( J)=0. 

HI Y  < J ) =0 i 

IF  < KOUT  EO:  I  OUT ) WRITE ( NWT  >  22 ) 

FORMAT ( /  /  r 10 X  f '  C  V  8  E  M  APPROXIMATION  FUNCTION  NODAL  VALUES • ' f // ? 
C6X» 'NODE ' ?6X» 'STATE' j 14X? 'STREAM' >/?3X» 'NUMBER' ?3X. 'VARIABLE' » 
C12X •  'FUNCTION' ) 

DO  1000  J= 1 r  NNOD 

CALCULATE  BOUNDARY  ELEMENT  CONTRIBUTIONS 


DO  500  K= 1 »  NNOD 
K  N  =  K  +  J. 

I F  <  KK  s GT . NNOD ) KKa  1 

IF'K.EQ. J.OR.KK.EQ. J)GOTO  500 

CALL  CAUCHK  J  •  K  »  KK  ?  A  ;  B  »  C  f  D  ) 

C  1  =  REX  < KK  >  *  (  X  (  J  >  -X <  K  )  >  -REY  ( KK  >  *  (  Y  '  J  >  -Y  <  K  )  ) 
C-REX(K)*(X( J)-X(KK>  )  +REY  <  K  )  ♦  (  Y  <  J  )  -  Y  (  KK  )  ) 
C2=REX<KK)*'Y(J)-Y(K)  )  +REY  <  KK  )  %  (  X  <  J)-X(K)  ) 
C -REX  <  K  )  ♦  (  Y ( J)-Y (KK) ) -REY ( K ) * < X ( J ) -X ( KK ) ) 
H1Y<  J'=H1X( J)+C1*C-C2*P 
H1Y(J'=H1Y(J>+C1*D+C2*C 
CONTINUE 

CALCULATE  PRINCIPLE  VALUE  CONTRIBUTIONS 
K  =  J  -  1 

IF(K ,LT : 1  )K  =  NNOD 
KK= J+ 1 

IF <  K  K  .  G  T . NNOD ) KK* 1 

YLN  =  SQRT( (Y(KK)-Y(J) >**2-MX(KK> - X  <  J ) )**2) 
XLM-SQRT ( (Y(K)-Y(J) )*#2+(X(K)-X(J) )**2) 
XXX*XLN/XLM 
A J=ALOG ( XXX ) 

AJ=DLOG(XLN/XLM) 

PJ= ( 360: -ANGLE ( J ) ) / 1 80 . *3 . 1 4 1 593 
H1X<  J)=H1X<  J)+REX<  J)*AJ-REY<  J)*BJ 
H1Y< J)=H1Y<  J)+REX( J)tBJ+REY< J)*AJ 

DIVIDE  BY  2fPI*i 


TEMP=H 1 X ( J ) 


non 


H 1 X  <  J  >  =  H 1  Y  <  J  '  /  6  . 29318 
HI Y ( J)*-TEMP/6 .29318 

IF ( KOUT  :EQ.  I0UT>URITE<NUT*450> JjH1X< J)*H1Y< J) 
450  F0RMAT<3X*I5*5X*F10,4f 10X*F10.4> 

1000  CONTINUE 


CALCULATE  NODAL  POINT  RELATIVE  ERROR 


2000 

200 


IF  <  KOUT  ,NE.  I OUT ) GO  TO  200 
URITE(NWT  *550) 

F0RMAT<///*10X» 'NODAL  POINT  RELATIVE  ERROR  VALUES''? 
DO  2000  1  =  1*  NNOD 
DA  =  REX  < I ) - H 1 X ( D 
0B=REY< I )-HlY(  I  ) 

WRITE(NWT  *450) I  *  DA » DB 
CONTINUE 

I F  <  KOUT  . EQ  >  999? ) KOUT  =  KP 

RETURN 

END 


H 
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SUBROUTINE  ANGLE 


SUBROUTINE  ANG<NNQD) 


CONHQN/BLK  l/XUOO) 
CONMON/BLK  2/YU00) 
COMMON/BLK  7/ANGLE' 100> 


THIS  SUBROUTINE  CALCULATES  THE  ANGLE  BETWEEN  EACH  NODAL  POINT 


IOO 


P I =ACOS ( - 1  .  ) 

DO  100  I  =  1 »  NNOD 
J=I-1 
J  J  =  I  + 1 

IF  < J.EQ.O) J  =  NNOD 
IF' JJ.GT.NNOD) JJ=1 
XJ=X( J)-X<I> 

XJJ  =  X( J  J  >  —X  < I > 

YJ=Y< J)-Y< I ) 

Y  J  J  =  Y ( J J ) - Y  ' I > 

CALL  C AUCH5 ' X JJ  *  Y  J  J  t  A  J J  > 

CALL  CAUCH5<XJ-YJf AJ> 

ANGLE'.  I  )  =  (  AJ-AJJ.' Y190  :  /PI 

IP  <  ANGLE  <  I )  .  L  T  :  0  •  )  ANGLE'  I  >=  ANGLE'.  I  >+360. 

CONTINUE 

RETURN 

END 
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SUBROUTINE  FRONT 


SUBROUTINE  FRONT ( NNOD t NFRS , NFRE r COND t XL AT , PQR , DELT t H 1 Y * KODB > 

THIS  SUBROUTINE  CALCUALTES  THE  NEW  INTERNAL  ANCLES  AND  NEW  POSITIONS 
OF  THE  FREEZING  OR  THAWING  FRONT  AFTER  EACH  TIME  INCREMENT  BY 
SHIFTING  THE  POSITIONS  VERTICALLY  OR  NORMALLY. 

COMMON/BLK  1/X<100> 

COMMON/ BLK  2/YC100) 

COMMON/BLK  7/ANGLE' 100) 

DIMENSION  0(50) rXP'50) .YP(50) 

DIMENSION  HIY(IOO) 

DO  50  1=1*50 
0  Q  < I ) =0 . 

...APPROXIMATE  T HE  EFFLUX  ALONG  THE  FREEZING  FRONT 
J  =  0 

DO  100  I=NFRS*NFRE-1 
J  =  J  +  1 

XX  =  X< 1  +  1  '  -  X  (  I  > 

YY  =  Y  < I  +  t )-Y< I) 

DIS  =  SQRT ( XXtXX  +  YYtYY  > 

FLUX=. 5*C0ND*(H1Y< I+1>-H1Y'I))/DIS 
O'  J  )  =  Q  J  >  +  F  L  U  X 
Q< J+l >=Q< J+l ) +FLUX 
00  CONTINUE 

UPDATE  THE  NEW  FREEZING  FRONT 
J  =  0 

DO  200  I =NFRS • NFRE 
J  =  J+1 

C ...  DETERMINE  THE  NORMAL  DIRECTION 
I P  1  =  I  +  1 

IF'IPl  : GT .  NNOD) IP  1  =  1 
I M 1  =  I  - 1 

IF ( I  Ml  >LT.  1 ) I M 1 =NNOD 
P I = ACOS ' - 1  :  ) 

XJ  =  X( IM1 )-X( I  ) 

Y J  =  Y ' IM1 > -Y  < I ) 

CALL  CAUCHS'XJ* YJ*AJ) 

C ...  CALCULATE  THE  NEW  FREEZING  FRONT 
DELS  =  G<  J)*DELT/(XLAT*POR) 

IF ( I . EQ . NFRS  .OR,  I , EQ . NFRE ) GO  TO  250 
ANGL=. 5* (360-ANGLE' I ) XPI/180 , +AJ 
IF ( KODE  .EQ.  2 ) GO  TO  220 
XP(I )*X( I) 

YP  < I >  =Y ( I ) -DELS 
GO  TO  200 

220  XP(I>=X(I)+DELS*COS(ANGL) 

YP( I >=Y( I )+DELS*SIN< ANGL) 

GO  TO  200 

250  YP(I)*Y'I)-DELS*2. 

XP(I)-X(I) 

200  CONTINUE 

DO  300  I«NFRS.NFRE 
X(I)»XP(I> 

Y(I)*YP(I) 
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